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Abstract 


This  document  describes  a  numerical  model  that  was  developed  to  study  two-dimensional, 
reduced-gravity,  shallow-water  flows.  When  the  dynamics  of  these  flows  is  strongly  non¬ 
linear,  the  flow  may  become  hydraulically  supercritical  and  discontinuities  in  the  flow  field 
may  arise.  The  presence  of  discontinuities  in  the  flow  field  requires  a  special  numerical  treat¬ 
ment  in  order  to  maintain  both  accuracy  and  stability  in  the  numerically-approximated  solu¬ 
tion.  In  this  model,  a  shock-capturing  scheme  called  the  Essentially  Non-Oscillatory  (ENO) 
scheme  is  implemented.  The  ENO  scheme  is  a  high-order,  adaptive-stencil,  finite-difference, 
characteristic-based  scheme  for  hyperbolic  equations  that  has  been  applied  widely  to  flows 
governed  by  the  Euler  equations  of  gas  dynamics.  The  model  described  in  this  document 
was  developed  for  geophysical  applications,  and  therefore  includes  the  effects  of  rotation 
(constant  Coriolis  parameter),  forcing  (time  dependent  and/or  spatially  varying),  and  bot¬ 
tom  drag  (linear  or  nonlinear).  The  presentation  includes  the  mathematical  formulation  of 
the  model  as  well  as  instructions  on  how  to  prepare  and  execute  model  runs. 
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1  Introduction 


The  numerical  model  described  in  this  document  was  developed  as  part  of  a  research  effort 
conducted  by  Roger  Samelson  and  myself  on  orographically-modified  winds  in  the  coastal 
marine  atmospheric  boundary  layer,  funded  by  the  ONR  Coastal  Meteorology  Accelerated 
Research  Initiative,  grant  N00014-93- 1-1369. 

The  model  was  developed  to  investigate  hydraulically  transcritical,  two-dimensional, 
reduced-gravity  shallow-water  flows  (Rogerson  1999) .  The  dynamics  of  these  flows  is  strongly 
nonlinear,  and  when  the  flow  velocity  exceeds  the  gravity-wave  phase  speed  the  flow  becomes 
hydraulically  supercritical  and  discontinuities  in  the  flow  field  may  arise.  The  presence  of 
discontinuities  in  the  flow  field  requires  a  special  numerical  treatment  in  order  to  main¬ 
tain  both  accuracy  and  stability  in  the  numerically-approximated  solution.  In  this  model, 
the  Essentially  Non-Oscillatory  (ENO)  shock-capturing  scheme  is  implemented.  The  ENO 
scheme  is  a  high-order,  adaptive-stencil,  finite-difference,  characteristic-based  scheme  for 
hyperbolic  equations  that  has  been  applied  widely  to  non-rotating  flows,  including,  for 
example,  flows  governed  by  the  2-D  Euler  equations  of  gas  dynamics.  For  the  present  appli¬ 
cation  to  the  reduced-gravity  shallow-water  system,  the  effects  of  forcing  (time  dependent 
and/or  spatially  varying),  bottom  drag  (linear  or  nonlinear),  and  rotation  (constant  Coriohs 
parameter)  have  been  included. 

The  model  is  designed  to  approximate  a  solution  on  a  user-specified  orthogonal  curvi¬ 
linear  grid.  The  original  application  involved  channel-Uke  domains  with  varying  side-wall 
geometries  which  could  be  discretized  onto  an  orthogonal  curvilinear  grid  (using  a  conformal 
mapping  program  developed  by  Wilkin  and  modified  successively  by  Signell,  by  Samelson, 
and  by  Rogerson).  In  this  case,  variations  in  the  side- wall  geometry  can  lead  to  hydraulic 
criticality  in  the  two-dimensional  flow.  In  contrast,  problems  in  classic  hydraulic  flow  theory 
typically  involve  a  rectilinear  channel  geometry  (or  one-dimensional  geometry)  with  varying 
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bottom  topography.  The  presence  of  bottom  topography  is  not  currently  included  in  the 
model,  but  could  be  incorporated  easily  to  simulate  flows  in  which  the  layer  interface  never 
intersects  the  bottom.  For  cases  in  which  the  layer  depth  goes  to  zero,  significant  modifi¬ 
cation  would  be  required  and  it  is  possible  that  another  approach  (i.e.,  using  a  different 
numerical  scheme  entirely)  would  be  more  fruitful. 

The  types  of  boundary  conditions  that  are  currently  implemented  in  the  model  refiect 
the  fact  that  a  channel  geometry  was  used  in  the  original  application,  although  it  is  also 
possible  to  specify  doubly-periodic  boundary  conditions.  For  the  channel  configuration, 
one  wall  of  the  channel  may  be  “opened.”  All  walled  boundaries  are  free-slip.  As  with 
all  numerical  models,  the  boundary  treatment  is  critical  to  the  stability  of  the  solution. 
Because  the  ENO  scheme  has  high  accuracy  and  low  numerical  dissipation,  it  is  particularly 
sensitive  to  inappropriate  and/or  inaccurate  boundary  treatments.  The  presence  of  rotation 
in  particular  requires  careful  consideration  in  terms  of  the  numerical  treatment  of  walled 
boundaries. 

This  document  has  been  created  with  the  hope  that  it  will  serve  as  a  useful  aid  to  re¬ 
searchers  who  want  to  use  and/or  modify  the  model.  The  presentation  is  fairly  technical,  in 
that  it  includes  a  complete  description  of  the  ENO  algorithm  and  some  of  the  mathematical 
formalism  behind  the  scheme,  as  well  as  instructions  on  how  to  prepare  and  execute  model 
runs.  The  model  equations  are  formulated  in  Section  2.  In  Section  3,  the  Essentially  Non- 
Oscillatory  (ENO)  scheme  is  introduced,  and  a  detailed  description  of  the  ENO  algorithm 
as  it  applies  to  the  flux  in  a  one-dimensional  scalar  equation  is  presented  as  an  example. 
The  application  of  the  ENO  algorithm  to  the  two-dimensional  shallow- water  system  is  dis¬ 
cussed  in  Section  4.  Finite-difference  approximations  for  the  non-conservative  terms  in  the 
model  equation  are  presented  in  Section  5,  followed  by  the  time-stepping  scheme  in  Sec¬ 
tion  6.  Boundary  conditions  and  their  implementations  are  presented  in  Section  7.  The 
steps  required  to  prepare  and  execute  the  model  are  outlined  in  Section  8.  An  auxiliary 
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program,  swgrid.f ,  that  can  be  used  to  generate  an  orthogonal  curvilinear  computational 
grid  is  described  in  the  Appendix. 
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2  Model  Equations 


Consider  the  rotating,  reduced-gravity  (1-1/2  layer),  shallow-water  flow  governed  by  the 
equations, 

ut  +  u^Vu-\r  f{k  xu)  =  -g^Vh - Vpa - r'^B  (1) 

p  pfi 

ht  +  V-{uh)  =  0  (2) 

where  u  =  (u,v)  is  the  horizontal  velocity  vector,  h  is  the  depth  of  the  layer  (which  is 
assumed  to  be  nonzero),  p  is  the  density  of  the  layer,  g'  =  g/^pjp  is  the  reduced  gravity,  / 
is  the  Coriolis  parameter,  pa{x,y,t)  is  the  imposed  pressure  forcing,  and  tb  is  the  stress  at 
the  bottom  of  the  layer.  The  bottom  stress  is  typically  parameterized,  and  in  the  present 
case  takes  the  form, 

Tb  =  pCd\u\u. 

Equations  (l)-(2)  are  nondimensionalized  using  length,  velocity,  time,  and  layer  depth 
scales  L*,U*,t*  =  L*/U*,  and  D*,  respectively,  to  yield, 

ut  +  u-Vu  +  fo{kxu)  =  -F~^'Vh-'VP -^\u\u  (3) 

ht  +  ^-{uh)  =  0  (4) 

where  /o  =  fL*/U*  is  the  inverse  Rossby  number,  =  g'D*/U*'^  =  is  the 

squared  inverse  of  the  scaling  Proude  number,  VP  =  {Vp/ p){L* /U*'^)  is  the  nondimen- 
sional  pressure  gradient  divided  by  the  density  of  the  layer,  and  r  =  CdL*/D*  is  the 
nondimensional  drag  coefficient. 

Equations  (3)-(4)  can  be  generalized  to  any  orthogonal  curvilinear  coordinate  system. 
If  (C,  t])  are  the  coordinates  in  the  orthogonal  curvilinear  system,  then  the  change  in  the 
position  vector  x  =  {x,  y)  in  the  Cartesian  system  can  be  written  as, 

dx  =  miS(^  ^  -I-  m^drj  t) 
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where  m\  and  m2  are  the  coordinate  metrics  given  by 


m,  =  ^(! 

(Batchelor  1967).  It  follows  that  the  gradient  of  a  scalar  <f>  is, 

4> 

yi,  7«2  V7I  j 

the  divergence  of  the  vector  v  —  {vi,V2)  is, 


^driJ  \dT]J 


\midC  m2  dr] 


V  • «  = 


1 


mi  m2  L  dC 


d{m2Vi)  ^  d{miV2) 


and  the  gradient  of  v  in  the  direction  of  n  is. 


n  •  Vv  =  C 
+  fi 


n  ■  Vvi  + 
n  •  Vv2  — 


V2 


mim2 

m\m2 


drj 


(  dm\  dm2  \ 

dm\  dm2  \ 

1X2 


(  dmi 


9c  jr 


(5) 

(6) 


Therefore,  in  the  (C,77)-coordinate  system  Equations  (3)-(4)  become. 


1  1  1  /  ^  ir 

Ut  H - uur  H - vUn  H - v{umi.  —  nm2,)  —  fov  = 

mi  ^  m2  mim2  ' 

1  1  1  /  ^  ,  X 

vt  H - uvr  H - Wf, - u{umi„  —  vm2j  +  fou  = 

mi  ^  m2  mim2  ’  ^ 

-±(F-\  +  Pr,)-l\u\v 

m2  ^ 

ht  H - —  [{m2uh)r  +  {mivh)j,]  =  0 

m\m2 

where  now  u  and  v  axe  the  velocity  components  in  the  ^  and  77  directions,  respectively. 


(7) 

(8) 
(9) 


Equations  (7)-(9)  can  be  rewritten  in  flux  form  in  terms  of  the  state  vector. 


f  uh\ 

fu\ 

vh 

= 

V 

\  h  J 

[  h  J 
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yielding 


where 


qt  +  i^[Fiq)]c  +  i[G(g)],  =  C  +  T  +  V  +  M 


f  U^lh  +  F-^h‘^/2 
F  =  UV/h 

I  U 


(  VVjh 
G  =  y  Vh  +  ^-2^2/2 

I 


are  the  fluxes  in  the  C  and  fj  directions,  respectively,  and 

f  foV  \  I  -hPjmi  \  (  -r\U\Ulh^ 

C  =  -foU  ,  T  =  -hPr,lm2  ,  V  =  -r\U\V/h^ 

\  0  I  \  0  )  \0 


M  = 


mi  m2 


-aiV/h-a2U/h 

aiU/h-a2V/h 


are  the  terms  resulting  from  the  Coriolis,  pressure  forcing,  bottom  stress,  and  grid-metric 
contributions,  respectively.  In  the  expression  for  Ad, 


=  U mi  —  V m2, ,  02  =  U m2,  -t-  V mi^ 
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3  The  Essentially  Non-Oscillatory  (ENO)  Scheme 


The  Essentially  Non-Oscillatory  (ENO)  scheme  is  a  numerical  method  for  hyperbolic  con¬ 
servation  laws  of  the  form, 


^  =0  (or  =  r (w,  x,t),  a  residual  or  forcing  term)  (11) 

2=1 

(12) 


dxi 

u{x,t  =  0)  =  Uo{x) 


where  u{x,  t)  =  (ui,  tt2, . . . ,  UmV ,  x  =  {xi,X2,...,  Xd)’^,  d  is  the  spatial  dimension  of  the 
problem,  m  is  the  number  of  independent  variables,  ajid  t  is  time.  The  system  is  hyperbolic 
in  the  sense  that  any  linear  combination  of  the  Jacobian  matrices, 

,  M 

*  du 

for  real  7i  =  (71,72,  •••  )7d))  always,  has  m  real  eigenvalues  and  a  complete  set  of  eigenvec¬ 
tors. 

It  is  well  known  that  the  solution  to  this  equation  can  develop  discontinuities  even  when 
the  initial  condition  wo(®)  is  smooth.  Traditional  finite-difference  techniques  will  provide 
poor  numerical  approximations  of  the  solution  in  this  case.  In  general,  shock-capturing 
schemes  aim  to; 

•  achieve  high  accuracy  in  regions  where  the  solution  is  smooth; 

•  Tnaintaiu  sharp  profiles  of  discontinuities  (i.e.,  avoid  excessive  numerical  dissipation); 

•  avoid  spurious  oscillations  in  the  vicinity  of  discontinuities; 

•  accurately  represent  the  location  and  speed  of  discontinuities;  and 

•  avoid  non-physical  solutions  (e.g.,  entropy- violating  expansion  shocks). 

The  ENO  scheme  satisfies  these  criteria.  In  addition,  the  ENO  scheme  is  globally  high-order, 
losing  only  one  order  of  accuracy  in  discontinuous  regions  compared  to  smooth  regions. 
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A  good  review  of  numerical  methods  for  conservation  laws  can  be  found  in  the  book 
by  LeVeque  (1990).  A  more  complete  description  of  the  ENO  scheme  can  be  found  in  the 
papers  by  Shu  and  Osher  (1988  and  1989),  and  references  contained  therein. 

3.1  ENO  Primer 

To  illustrate  the  fundamental  principles  of  the  ENO  scheme  and  simplify  the  discussion, 
consider  the  one-dimensional  scalar  equation, 

^t  +  [/(^)]i  =  0 

for  some  function  f{u). 

A  numerical  approximation  to  (13)  is  called  conservative  if  it  is  of  the  form, 

-  /.-I/2)  (14) 

where  h+\/2  =  /(ui-i, •  •  •  ,«i+fc)  for  some  l,k  >  0  (LeVeque  1990).  The  key  to  shock¬ 
capturing  schemes  hinges  on  how  the  numerical  fluxes  at  the  “half”  grid  points,  fi+i/2j 
estimated.  By  using  a  conservative  method,  the  Lax-Wendroff  theorem  guarantees  that  if 
the  numerical  scheme  converges  to  a  numerical  solution,  then  the  numerical  solution  does 
indeed  approximate  a  weak  solution  of  the  partial  differential  equation  (Lax  and  Wendroff 
1960).  Convergence  of  the  numerical  scheme  can  often  be  proved  if  the  scheme  is  total- 
variation  diminishiTig  (TVD)  or  total- variation  bounded  (TVB)  (Harten  1984).  The  total 
variation  is  defined  as, 

TV{u)  =  ^  |«i+i  -Ui\ 

i 

and  the  scheme  is  termed  TVD  if 

TV(u”+^)  <  rv(u") 

for  all  n.  The  scheme  is  termed  TVB  in  0  <  t  <  T  if 


TV(u”)  <  B 


for  some  fixed  B  that  depends  only  on  TV  and  for  all  n  and  At  such  that  0  <  nAt  <  T. 
The  ENO  scheme  is  a  descendant  of  TVD  and  TVB  schemes,  but  is  unique  in  that  it  uses 
adaptive  stencils  to  compute  the  numerical  approximations  of  the  fiux,  fi^i/2  ■  High  accuracy 
is  achieved  by  approximating  /i+1/2  and  /i_i/2  to  very  high  order  in  such  a  way  that  the 
difference  yields  a  high-order  estimate  for  the  derivative  df/dx  (as  opposed  to  seeking  a 
high-order  approximation  for  dfjdx  directly).  The  use  of  an  (r  +  l)-point  adaptive  stencil 
yields  (r  -I-  l)-order  accuracy  in  smooth  regions  of  the  flow  and  r-order  accuracy  right  up  to 
discontinuities,  and  is  formally  r-order  accurate.  A  high-order  interpolating  polynomial  is 
constructed  at  each  time  step  to  approximate  the  flux  from  information  at  the  surrounding 
grid  points,  avoiding  regions  where  discontinuities  are  present.  For  example,  to  compute 
fi+i/2  using  a  4-point  stencil  in  a  smooth  region  of  the  flow,  the  centered  stencil 

^1+1/2  .f  (ui— 1,  Uj,  ^2.^1 ,  Uj.1-2) 

would  be  utilized,  while  in  the  presence  of  a  local  discontinuity,  say  located  near  Xi-i,  the 
stencil  would  be  shifted  to  the  right  to  obtain  information  from  the  smoother  region,  e.g., 

fi+l/2  =  /(Ui)Ui+i,U2+2,U2+3). 

The  fact  that  the  scheme  involves  an  adaptive  stencil  application  has  hindered  progress 
towards  a  convergence  theory  for  the  ENO  scheme.  Nevertheless,  numerical  convergence 
has  been  demonstrated  empirically  in  a  number  of  challenging  problems  in  gas  dynamics 
including  Riemann  problems,  shock-wave  interactions,  and  shock-turbulence  interactions 
(see,  for  example,  Hannappel,  Hauser  and  Friedrich  1995). 

To  aid  the  stencil  selection  process  and  the  construction  of  the  interpolating  polynomial, 
divided  differences  are  computed  for  the  flux.  The  divided  differences  of  a  function  w  =  w{x) 
are  defined  recursively  as, 

w[Xi]  =  w{Xi) 
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w[xi, .  .  .  ,Xi+k] 


where  w[xi, . . . ,  Xi+k]  denotes  the  divided  difference  of  w  of  order  k.  If  w  is  smooth  (i.e.,  w 
is  infinitely  differentiable;  w  €  C°°)  in  the  interval  [xi,Xj+fc],  then 


1  ^i+/c] 


1 

A:!  dx’^  ’ 


but  if  w  is  discontinuous  in  the  p-th  derivative  (0  <  p  <  A;),  then 


w[xi, Xi+k]  =  O  (As 


where  denotes  the  jump  in  thep-th  derivative  (Harten  et  al.  1987).  Therefore,  divided 
differences  can  be  used  to  detect  discontinuities.  Now  consider  the  function  h{x)  such  that, 

1  rx-^Ax/2 

/(u(x))  =  y-  l  ,  HO<K. 

/\X  Jx—Ax/2 

It  follows  easily  that  if  H  is  the  primitive  function  of  h,  i.e., 

H{x)  =  r  hiOdC 
J—OO 

then 

H{xi+ij2)=  r^^'^  h{CjdC=  ^  h(C)dC  =  ^  (a:fc+i/2-a:fc_i/2)/K) 

•^-°o  ,  fc=-oo''^fc-V2  k=-oo 

and 

H{xi^il2)  -  i?(xj_i/2)  =  (a:i+i/2  -  ^i-il2)f{ui)- 
Therefore,  the  divided  differences  of  H  can  be  obtained  directly  from  the  divided  differences 
of/, 


H[Xi-ll2-,Xi-ll2+\\  —  fWi] 
H[xi^i/2,  •  •  • ,  a:i-i/2+fe]  =  ■  ■  ■  ’ 
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3.1.1  Basic  ENO-Roe  Algorithm 


We  seek  a  solution  to 

+  [/  (^)]x  —  0 

in  which  [/(tt)]x  is  approximated  via, 

^  (/i+i/2  -  /i-1/2)  • 

The  numerical  fluxes  /  are  computed  to  r-th  order  using  the  algorithm  outlined  below.  To 
simplify  the  notation,  we  will  denote  the  A:-th  divided  difference  of  H  at 


^i—l/2  —  -^^[^1—1/21  ■  •  ■  1/2+A:]- 


The  “ENO-Roe”  algorithm  (Shu  and  Osher  1989)  for  /i+1/2  is: 


1.  Compute  the  divided  differences. 


Hi-U2 

^i-l/2 


=  /N  =  f{ui) 

~  •  •  •  )  1]»  fc  =  2, . . .  ,  7" -t- 1 


2.  Estimate  the  local  sign  of  df/du  at  by  computing  the  Roe  speed  (Roe  1981), 

fiui+i)  -  f{Ui) 


Oi+l/2  — 


Ui+l  Ui 


3.  Let  s{k)  denote  the  starting  stencil  point  at  stage  k  in  the  selection  process.  Select 
the  first  stencil  point  {k  =  1)  in  the  upwind  direction, 


s(l) 


_{  i  if 
1  i  +  1  otherwi 


+1/2  >  0 


otherwise 

4.  At  each  stage  thereafter,  add  an  additional  point  to  the  stencil  from  the  “smoother” 
region,  using  the  difference  table  for  the  comparison; 

s{k)  —  1  if  Rs(fc)_i/2  ^  ■^s(fc)-l/2-l 


s{k  +  l) 


=1 


s{k)  otherwise 
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When  starting  stencil  location  is  backed  up;  we  add  the 

point  to  the  “left”  of  previous  starting  point.  When  H^(^k)-i/2  < 
starting  stencil  location  is  unchanged;  we  effectively  add  a  point  to  the  “right-hand” 
end  of  the  stencil. 

5.  After  the  (r-f-l)-point  stencil  has  been  selected,  we  construct  a  high-order  interpolating 
polynomial  for  the  primitive  function, 


+  •^5(2)-1/2(®  ~  ®5(1)-1/2)(^  ~  ^s(l)-l/2-i-l) 

-h  ■ffs(3)_i/2(a;  -  Xs^2)-1/2){x  -  Xs{2)-l/2+l){^  -  a:s(2)-l/2+2)  +  •  •  • 


r+l 


rk-i 


""  ^s{l)-l/2)  +  I  ^s{k)-l/2 


n(^  -  a:s(fc_l)-l/2-f-a) 

-Q:=0 


whose  derivative  is, 


+  H'^{2)-i/2[i^  -  ^4(i)-i/2)  +  -  a:s(i)-i/2-t-i)]  +  •  •  • 


rH-1 


=  +  E 

k=2  [ 

The  interpolating  polynomial  for  /i+1/2  is  then, 


■fc-i  fc-i 

E  n  {x  -  Xs{k-l)-l/2+p) 

q;=0  0=0, a^P 


/.«/2  = 


■i-fl/2 


dx 

=  ■^s\l)-l/2 

r+1  f 

+  X)  S  Hs{k)-l/2 

k=2  \ 


k-l  ifc-1 

E  n  (^i+1/2  -  ^s(k-l)-l/2+p) 

[_a=0  0=O,a-^0 


(15) 


3.1.2  Entropy  Fix;  the  ENO-LLF  Algorithm 

The  ENO-Roe  scheme  described  above  does  admit  a  non-physical  entropy-violating  expan¬ 
sion  shock  but  the  problem  can  be  easily  remedied  (Shu  and  Osher  1989).  If  f'[u)  does 


12 


not  change  sign  between  Ui  and  tti+i,  then  we  compute  the  numerical  flux,  fi+1/2,  accord¬ 
ing  to  Equation  (15)  in  the  ENO-Roe  fashion.  If  f'{u)  does  change  sign  between  Ui  and 
Uj+i,  then  we  compute  the  numerical  flux  in  a  slightly  difierent  fashion  based  on  the  local 
Lax-Priedrichs  flux  (the  “ENO-LLF”  scheme)  described  below. 

The  flux,  /(«),  can  be  split  into  two  parts 


where 


with 


f{u)=f^{u)  +  f  {u) 

/■^(“)  =  \[fiu)  +  au] 

/“(«)  =  \[fiu)-au] 

a  =  max  |/^(u)| 


so  that 


df+ 

du 

du 


>  0 


<  0. 


The  numerical  flux  is  similarly  split  in  the  Lax-Friedrichs  fashion. 


fi+1/2  —  /i+1/2  fi+1/2 

and  ENO  approximations  axe  computed  for  each  component  as  follows; 

1.  Compute  the  divided  diflerences  for  '^H  and  ~H, 

^Hl_i/2  =  \{f[ui]±Oli^l/2u[Xi]) 

^^i—1/2  ~  2  ■  ■  ■  ’  ^  ■  ■  ■  ’  fc  =  2, . . . ,  r -f- 1 
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2.  Define 


Q!i+i/2  =  |/'(u)| 

3.  Select  the  first  stencil  point  for  the  +  and  -  components  in  the  upwind  direction  with 
respect  to  the  half-grid  point  i  -|- 1/2.  Since 

a/+i 


du 

dr 


^i+l/2 


du 

l^t+l/2 

the  first  stencil  points  are  chosen  as, 


—  5  Cii+l/2  >  0 


=  -\  0!j+i/2  <  0 


s+(l)  =  i 
s“(l)  =  i -I- 1 


4.  Select  the  rest  of  the  stencil  in  the  ENO  fashion, 


s^{k  +  l)  =  <  + 


3±(i)  -  1  if  >* 


$^{k) 


otherwise 


5.  Form  the  interpolating  polynomials  for  fi+i/2 


ft-.,.  = 


'i+l/2 


±  trl 

s±{l)-l/2 


r+1  ( 

+  Er-^?*w-v2 


k-1  fc-1 

E  n  {^ii-l/2  -  3^s±(fc-l)-l/2+^) 

7=0/?=0,7#/3 


(16) 


6.  The  numerical  fiux  computed  using  the  ENO-LLF  scheme  is  then. 


fi+l/2  —  f^l/2  -^1+1/2 


(17) 


To  prevent  entropy- violating  expansion  shocks  in  the  EN 0-Roe  scheme,  the  numerical 
fiux  must  be  computed  according  to  Equations  (16)-(17)  if  /'(«)  changes  sign  between  Ui 
and  Ui+i,  and  not  in  accordance  with  Equation  (15). 


Since  the  entropy  fix  for  the  ENO-Roe  scheme  requires  the  implementation  of  an  addi¬ 
tional  scheme,  ENO-LLF,  one  might  wonder  if  it  would  be  better  to  compute  all  of  the  fiuxes 
using  the  ENO-LLF  scheme  in  the  first  place.  Employing  the  ENO-LLF  scheme  globally 
would  certainly  be  simpler  (algorithmically)  than  ENO-Roe  with  entropy  fix,  and  in  the 
shallow-water  model  the  user  has  the  option  to  select  either  scheme.  However,  the  numerical 
dissipation  associated  with  the  ENO-Roe  scheme  is  less  than  that  of  ENO-LLF  (Shu  and 
Osher  1989),  so  there  is  less  shock  smearing  and  better  overall  accuracy  with  ENO-Roe. 
In  general,  I  have  found  the  ENO-Roe  (with  entropy  fix)  solutions  to  be  superior  to  those 
generated  by  ENO-LLF. 

3.1.3  Hybrid  ENO;  Biased  Stencil  Selection 

Adaptive  stencil  selection  is  the  key  feature  of  the  ENO  scheme.  It  allows  an  interpolating 
polynomial  to  be  constructed  using  a  stencil  that  avoids  discontinuous  regions  of  the  flow. 
It  is  inevitable  that  in  this  process  linearly  unstable  stencils  will  be  selected.  In  general, 
however,  the  selection  of  hnearly  unstable  stencils  does  not  lead  to  nmnerical  instability 
since  rapid  stencil  switching  is  often  observed  (Harten  et  al.  1987).  However,  in  smooth 
regions  of  the  flow,  the  use  of  linearly  unstable  stencils  (and  the  rapid  stencil  switching  that 
accompanies  it)  can  degrade  the  convergence  rate  of  the  solution  (Rogerson  and  Meiburg 
1990).  The  error  reduction  diuing  mesh  refinement  is  not  uniform  and  in  some  cases  grid 
refinement  can  produce  an  increase  in  the  truncation  error.  This  degeneration  in  acciuacy 
can  be  remedied  by  using  fixed  linearly  stable  stencils  in  smooth  regions  of  the  solution  and 
adaptive  stencils  where  strong  gradients  are  present.  A  simple  modification  to  the  basic 
ENO  algorithm  combines  the  use  of  fixed  and  adaptive  stencils,  creating  a  “hybrid”  ENO 
scheme  that  restores  the  desired  accuracy  (Shu  1990).  In  the  stencil  selection  process  (item 
4  in  ENO-Roe  and  ENO-LLF  algorithms)  we  simply  replace, 

s(k  -I-  1)  -  I  ^  ^s(A:)-l/2  ^  ^s{k)-l/2-l 

s[K  -h  i;  -  <  Otherwise 
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with 


if  s{k)  >  c{k)  then 


s{k  +  l)  = 


s{k)  —  1  if  ^  ^s{k)-l/2-l 

s{k)  otherwise 


else 


(18) 


s{k  + 1) 


_  J  -  I  if  -^s(fc)-l/2  —  ^^s(k)-l/2-l 

^  s{k)  otherwise 

where  c{k)  is  the  leftmost  grid  point  in  the  centered  stencil.  The  weighting  factor  of  2 
in  Equation  (18)  is  used  for  the  reasons  provided  by  Shu  (1990).  Restated,  the  modified 
algorithm  is, 

if  the  stencil  is  to  the  right  of  the  centered  stencil  [i.e.,  s{k)  >  c(A;)]  then 
favor  adding  a  point  on  the  left  [i.e.,  s{k  +  1)  =  s{k)  —  1] 
else 

favor  adding  a  point  on  the  right  [i.e.,  s{k  +  1)  =  ^(A:)]. 

The  modified  algorithm  biases  the  stencil  selection  towards  the  linearly  stable  centered 
stencil  in  smooth  regions  where  ^.nd  H’^(k)-il2-i  same  order  of  magnitude. 

3.1.4  Implementation  Issues 

Notice  that  on  an  equally-spaced  grid.  Equation  (15)  becomes, 


r+l 


/i+1/2  =  H]^l)-\j2  +  51  1  -^i(fc)-l/2(^^) 


k=2 


5.  k _ 1 

''■"E  n  (i-s(A:-l)-H-;0) 

a=0  /5=0,a^^ 


Therefore,  if  we  compute  the  undivided  differences, 


.  (19) 


^i-l/2  ~  /(‘^i) 


H: 


i-l/2 


=  n: 


k-l 


i+1/2  ^i-1/2 


yk-1 


A:  =  2, . . . ,  r  -I- 1 


(19)  becomes, 

r+l  ( 

/i+1/2  =  ^l(l)-l/2  +  E  1  '^i(A:)-l/2 

k=2 

Since  the  coefficients  in  the  summation. 


k-l  k-l 

^  n  {i-sik-i)+i-p) 

q:=0 


k-l  k-l 

E  n  (*  ”  ■®(^  “  f)  +  f  - /^)’  A:  =  2, . . .  ,r -I- 1 

a=:0  /3=O,a^j0 


(20) 
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depend  on  (i  —  s{k))  (the  difference  between  the  grid  point  in  question  and  the  left-most 
stencil  point)  and  not  on  i  itself,  the  set  of  possible  coefficients  can  be  precomputed  for  use 
in  Equation  (20).  A  similar  simplification  can  be  applied  to  the  split  numerical  fluxes  in 
the  ENO-LLF  algorithm  (Equation  (16)). 

When  the  grid  is  not  uniform,  we  mahe  a  change  of  variables,  e.g.,  x  -t  C?  reformulate 
the  governing  equation, 

“t  +  ^[/(“)]c  =  0 

m 

where  m  =  dx/dC  is  the  grid  metric.  The  approximation  for  the  flux,  fi+1/2,  then  proceeds 
on  the  C  grid  for  which  AC  is  constant. 

3.1.5  Implemented  ENO  Algorithm 

The  ENO  algorithm  that  is  implemented  in  the  model  is  a  hybrid  EN 0-Roe  scheme  with 
the  “entropy  fix”  (Shu  and  Osher  1988;  Shu  and  Osher  1989;  Rogerson  and  Meiburg  1990; 
Shu  1990).  Below,  we  recapitulate  the  algorithm  for  clarity  as  it  applies  to  the  1-D  scalar 
equation  (13). 


1.  Compute  the  undivided  difierences  for  /  and  it. 


«Li/2 

njk 

^i-1/2 


=  f[Ui]=f{Ui) 
=  n[xi]  =  u{xi) 


^/[iti, ...  I  Uj-i-fc—l] 
—u\Xi,  .  .  .  ,  Ij+fc—l], 


fc  =  2, . . . ,  r  -h  1 


2.  Compute  the  Roe  speed. 


and 


<*i-t-l/2  — 


fjuj+i)  -fjuj) 


0‘i+i/2=  !/'(«)! 
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3.  Form  the  undivided  diflference  tables  for  the  split  local  Lax-Priedrichs  flux, 

^  «i+l/2  ^-l/2)>  fc  =  1, . .  -  , r  +  1 


4.  If  f'{u)  does  not  change  sign  between  Ui  and  Ui+i,  then  construct  /i+i/2  using  the 
ENO-Roe  algorithm. 


Select  the  first  stencil  point  in  the  upwind  direction, 


^  M  if  ai+i/2  >  0 
^  '  I  i  +  1  otherwise 


Select  the  remaining  stencil  points.  Bias  the  stencil  selection  towards  the  linearly 
stable  centered  stencil,  c{k),  in  smooth  regions, 
if  s{k)  >  c{k)  then 

s(k  + 1)  =  I  ^ 

'  1  s(A:)  otherwise 


s{k  + 1) 


^  r  S{k)  -  1  if  >  2^"(fc)-l/2-l 

I  s{k)  otherwise 


Compute  the  interpolating  polynomial  for  /i+i/2) 


/ii+l/2  =  ^s(l)-l/2 

r+1  C  FA:-!  fc-1 

+  ^  s  ^s(fc)-i/2  n  (i  —  s{k  -  1)  +  1  -  0) 

k=2  I  a=0  /3=0, 05^/3 


5.  If  f{u)  changes  sign  between  Ui  and  Ui+i,  then  construct  /i+1/2  using  the  ENO-LLF 
algorithm. 


Select  the  flrst  stencil  point  in  the  upwind  direction, 


s+(l)  =  i 


s-(l)  =  i  +  1 
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•  Select  the  remaining  stencil  points,  s'^{k)  and  in  the  hybrid  ENO  fashion, 


if  s^{k)  >  c{k)  then 

+  I)  =  f  »*('=)  -  1  if  2(*w»V)-i/2)  a* 
[  s^{k)  otherwise 

else 

s*(k  +  11  -  !  ^  “  Xt)-l/2  a 

+  otherwise 

•  Form  the  interpolating  polynomials  for  f^i/2  /i+i/2’ 


f±  — 

h+1/2  ■“  ^s=t(l)-l/2 

r+l  ( 

+  Er’^‘*w-V2 

fc=2  [ 


E  n  {i-s^{k-i)+i-p) 

a=0  ^=0,a^p 


Sum  the  split  fluxes  to  obtain  the  interpolating  polynomial  for  the  numerical 
flux, 

/t+l/2  =  /i+l/2  +  fi+1/2 


6.  Repeat  the  previous  steps  to  compute  fi-1/2  and  approximate  {df/du)i  as. 
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du 


3.2  Time-stepping  Scheme 


ENO  spatial  approximations  are  typically  paired  with  TVD  time-stepping  schemes  (see  Shu 
and  Osher  (1988)  for  background).  Shu  and  Osher  (1988)  formulated  several  Runge-Kutta 
time-stepping  schemes  that  have  optimal  (i.e.,  large)  CFL  restrictions.  For  the  equation, 

du 

the  2nd-order  and  Srd-order  Runge-Kutta  methods  are, 


u°-  =  u"  -1-  AtL{u^) 

(21) 

“  iir”  +  if- + 

(22) 
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and 


(23) 

(24) 

(25) 


=  u^  +  AtL{u^) 

yb  =  L"  +  1  a  1 
4  4  4 

The  theoretical  CFL  coefficient  for  both  schemes  is  1.  In  practice,  the  recommended  maxi¬ 
mal  CFL  coefficient  is  0.6  when  L{u)  is  approximated  with  an  ENO  algorithm  (C.-W.  Shu, 
personal  communication),  i.e., 

max  |/'(u)|  <  0.6.  (26) 

At.  V,  - 
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4  ENO  Scheme  for  the  Shallow  Water  System 

Reconsider  the  shallow- water  system  in  flux  form  on  a  curvilinear  grid  (see  Equation  (10), 
Section  2), 

Qt  +  ;;“[-^(9)]c  +  =  Q  (27) 

m\  m2 

where  ,  ,  ,  , 


q=z  \  vh  =  \  V 


U^/h  +  F-^h?/2 
UVIh 
U 


f  UV/h 

G  =  FV/i  +  F-^h'^/2 

I  F 


and  Q  =  C+F+V+Ai  is  the  sum  of  all  of  the  terms  on  the  right-hand  side  of  Equation  (10). 

For  systems  of  equations,  the  ENO  algorithm  is  applied  to  each  local  characteristic  field, 
not  to  each  state  variable.  To  illustrate  how  the  fluxes  in  the  ^  direction  are  computed, 
consider  the  one-dimensional  conservation  equation. 


9t  +  — [■F’(q)]c  =  0 
mi 


where 


Qt  +  —Aq^  =  0 
mi 


„„  f  2u  0  —u^  +  F-^h 
^  oF 

A  =  -;r-  =  V  U  —UV 

^9.10  n 


The  matrix  A  has  eigenvalues  and  left  and  right  eigenvectors. 


=  u 

1(1)  = 

( 

0  1 

-V  ) 

r(i) 

=  (01 

0 

)T 

A(^)  =  u  +  c 

1(2)  = 

-Tc  ( 

-1  0 

u  —  c  ) 

r(2) 

=  (  +  C  V 

1 

r 

1 

II 

1(3)  = 

h  ( 

-1  0 

u  +  c  ) 

r(3) 

=  (  —  C  V 

1 

r 

where  c  =  y/Fr  ^h.  Equation  (29)  can  be  projected  onto  the  eigenspace  via 

S-^qt  +  —{S-'^AS)S-^q^  =  0 
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where 


yielding, 


where 


s=( 


r»(l)  7*(2)  7*(3) 


). 


f  \ 

5-1  =  Z(2) 

V 


Rt  H - A.Rf 

mi 


R  = 

(  ""  ^ 

It  -|-  2c  1  , 

A  =  5-1 A5  = 

{  Ad) 

A(2) 

\ 

^  u  —  2c  j 

V 

Ad)  ) 

As  in  the  1-D  scalar  case  (see  Equation  (14),  Section  3.1),  the  spatial  derivative  in 
Equation  (28)  is  computed  using  the  simple  difference  formula, 
dF{q)\ 


dC 


^  {F i+l/2  -  Fi-1/2)  _  p 

(Ci+1/2  —  Ct-1/2) 


i+1/2  “  F i-if2 


where  .Fj+1/2  and  are  high-order  approximations  to  the  flux  obtained  from  an  adap¬ 

tive  stencil  that  avoids  discontinuous  regions  of  the  flow.  To  compute  Fi+1/2,  the  algorithm 
outlined  in  Section  3.1.5  for  the  1-D  scalar  case  is  generalized.  First,  undivided  difference 
tables  are  computed  for  each  component  of  the  flux  and  the  state  vector,  as  in  step  1  in 
Section  3.1.5,  i.e.. 


'^i-l/2  =  F[q,] 

1 

to 

II 

^tl/2  = 

^i-1/2  ~  ■  ■ 

The  difference  tables  are  then  projected  onto  the  eigenspace  using  the  left  eigenvectors  of 
A, 

Ri-\/2  =  ^)i+l/2'R'i-l/2- 


Only  the  portion  of  the  difference  table  that  might  be  utihzed  in  the  approximation  for 
-^i-t-1/2  is  projected  onto  the  eigenspace  (see  Figure  1).  For  each  projected  component,  p. 
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i— 3  i—2  i — 1  i  i+1  i+2  i+3  i+4 

i+]  /2 


Figure  1:  The  portion  of  the  diflFerence  table  that  could  possibly  be  utilized  in  the  approx¬ 
imation  of 
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an  appropriate  stencil  is  selected  in  the  ENO  fashion  and  the  interpolating  polynomial, 
formed,  as  outlined  in  steps  2-5  in  Section  3.1.5.  The  interpolating  polynomials 
in  the  eigenspace  are  then  projected  back  using  the  right  eigenvectors  of  A,  i.e., 

^i+ll2  =  Si+i/2F  i+i/2- 

The  Roe  speeds  in  this  case  axe  the  eigenvalues, 

„(p)  _  \{p)  „  _  1  0  Q 

and  the  local  Lax-Priedrichs  estimate  is. 


2+1/2 


ar.'wo  =  max 


p  =  l,2,3. 


Qi<q<Qi+i 

Since  the  Roe  speeds  and  projection  matrices  are  evaluated  at  the  half-grid  points,  a  suitable 
average  must  be  computed  for  u,  v,  and  h.  The  appropriate  averages  are  the  Roe-averaged 
quantities  (Roe  1981)  given  by. 


Ui+l/2  = 
‘>Ji+i/2  = 

hi+i/2  = 


Ujy/hj  -|- 

Vhi  +  \/hi+i 

Viy/Tn  +  Vj+ly/hj+l 

+  V^i+l 

-{hi  +  hi+i). 


The  ENO  scheme  is  easily  generalized  to  multi-dimensions  since  the  approximations  to 
the  fluxes  F{q)  and  G{q)  are  computed  separately.  Equation  (27)  is  therefore  approximated 
as 

Qt  =  -^(■^i+l/2J  -  F i-i/2,j)  -  —{Gij+i/2  -  Gij-1/2)  +  Qij-  (30) 

The  computation  of  (jjj+1/2  is  analogous  to  that  of  Em  completeness,  it  is 

noted  that  the  matrix 

V  u  —uv 

0  2v  -v^  +  F-'^h 
0  1  0  / 
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has  eigenvalues  and  left  and  right  eigenvectors, 


A(i) 

=  V 

= 

( 

1 

0  -u  ) 

r(i) 

=  ( 

1 

0  0 

A{2) 

=  V  +  C 

Z(2)  = 

1 

2c 

( 

0 

—  1  V  —  c  ) 

r(2) 

=  ( 

u 

V  +  c  1 

A(3) 

=  V  —  C 

Z(3)  = 

1 

2c 

( 

0 

—  1  V  +  c  ) 

r(3) 

=  ( 

u 

V  —  c  1 

)T 
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5  The  “Right-hand  Side” 


This  section  describes  the  finite-difference  approximations  for  the  non-conservative  terms 


on  the  right-hand  side  of  Equation  (10), 


q,  +  =  C  +  V  +  V  +  M 


where 


uh\  U\ 
q=  vh  \=  V  \ 

\  h  )  \h  J 

and  C,V,T>,  and  M  represent  the  Coriolis,  pressure  forcing,  bottom  stress,  and  grid-metric 
contributions,  respectively.  The  approximations  for  C,  P,  X>,  and  Ad  are  all  straightfor¬ 
ward. 

5.1  Coriolis  Term 


The  approximation  for  the  CorioUs  term  is  simply, 

f  foVij 
(^ij  ~  foUij 
\  0 


5.2  Pressure  Forcing  Term 


The  pressure  forcing  term  is, 


/  —hP(^/mi 
V  =  -hPrj/m2 

I  0 


The  model  is  currently  set  up  to  read  from  an  input  file  the  steady  component  of  the 
(nondimensional)  pressmre-gradient  forcing, 

^  dPHCv)  1  dP^i<:,v)\ 


VP^  = 


mi  dC  ’  7772  dr] 


along  with  the  initial  flow  fields.  If  there  is  a  time-dependent  component  to  the  pressure 


gradient. 


VP'  = 


1  ap'(c,r?,t)  1  dP'{c,v,t) 

mi  dC  ’  m2  dr] 
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it  should  be  specified  in  subroutine  Forcing.  (In  the  current  version,  VP'  corresponds 
to  a  propagating  low-pressure  anomaly.)  The  approximation  for  the  pressure  term  is  then 


simply, 


"P  ij  — 


/  -hiP^+Pl)/mi  \ 
-h  (pO+P')/m2 
0 


^  ij 


5.3  Bottom  Stress  Term 
The  bottom  stress  term  is, 

(  -r\U\U/h^  \ 

V  =  (  -r\U\V/h‘^  1 

where  r  =  CdL*ID*  is  the  nondimensional  drag  coefiicient.  The  appropriate  formula  for 
Cd  is,  of  course,  application  specific.  For  the  original  application  to  the  marine  atmospheric 
boundary  layer  flow,  for  example,  two  formulas  for  Cd  for  the  air-sea  interface  were  coded. 
The  first  is  the  10-m  neutral  drag  coeflBcient  given  by  Large  and  Pond  (1981), 


IO^Cd 


=li 


14  uio  <  10  m  s  ^ 

49  +  0.065uio  uio  >  10  m  s”’^. 


The  second  is  the  drag  coefficient  obtained  from  the  Coastal  Waves  96  field  experiment 
(Edwards  and  Rogerson,  in  preparation). 


IO^Cd 


={«: 


43  —  0.261uio  uio  <  6.1  m  s  ^ 

44  +  0.065uio  uio  >  6.1  m  s“^. 


In  both  cases,  the  10-m  winds  in  the  formula,  uio,  must  be  related  to  the  model’s  nondi¬ 
mensional  layer-averaged  wind  speed,  lu|.  This  relation  is  currently  specified  as. 


uio  =  0.75|ti|i7*. 


Note  that,  to  obtain  the  nondimensional  drag  coefficient,  r,  the  length,  velocity,  and  depth 
scales  {L*,  U*,  and  D*)  must  be  specified  as  well. 
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A  linear  bottom  stress  may  also  be  modeled.  In  this  case, 


/  -rU/h 
-rV/h 


\  0 


where  f  is  a  (constant)  nondimensional  linear  drag  coefficient  provided  as  model  input. 
5.4  Grid  Metric  Term 


The  grid  metric  term  is, 


M  = 


1 

mim2 


-aiV/h  -  a2U/h  \ 
aiU/h-a2V/h 
-012 


where 


ai  =  Um\^  —  Vm2^,  02  =  Um2^  + 


The  metric  derivatives. 


dmi  dm2 


are  computed  during  the  initialization  phase  of  the  model  using  simple  first-order  differences. 
The  approximation  for  Mij  during  the  model  integration  follows  directly  from  the  algebraic 
expression  above. 
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6  Time-stepping  Scheme 


The  numerical  solution  to 

Qt  =  L{q)  =  -;^[^’(g)]c  -  +  S  (31) 

is  advanced  in  time  with  a  TVD  Runge-Kutta  scheme  (see  Section  3.1  for  the  definition  of 
TVD).  In  the  shallow- water  model,  the  user  can  specify  either  2nd-order  or  3rd-order  TVD 
Runge-Kutta  time-stepping  in  the  form, 

=  g”-hAiL(g”)  (32) 

qn+r  =  ig"  +  ig“-hiAtL(g“)  (33) 

or 

=  g"-hAtL(g”)  (34) 

=  |9”  +  |9»  +  iA(L(9“)  (35) 

«”+'  =  +  +  (36) 

respectively.  As  mentioned  in  Section  3.2,  the  theoretical  CFL  coefficient  for  both  schemes 
is  1.  In  practice,  however,  the  recommended  maximal  CFL  coefficient  is  0.6,  i.e., 

when  L{u)  is  approximated  with  an  ENO  algorithm  (C.-W.  Shu,  personal  communication). 
The  use  of  lower  CFL  coefficients  (e.g.,  0.1  or  0.2)  is  frequently  quoted  in  the  literature  as 
well.  For  the  shallow-water  model,  I  have  typically  selected  CFL  coefficients  in  the  range 
0.4-0.5  and  would  categorize  the  use  of  CFL  coefficients  in  the  range  0. 1-0.2  as  conservative. 


7  Boundary  Conditions 

The  types  of  boundary  conditions  that  are  currently  implemented  in  the  model  correspond 
to  two  basic  geometries:  a  channel-like  domain  and  a  doubly-periodic  domain.  For  the 
channel  configuration,  the  along-channel  direction  is  assumed  to  be  the  77  direction.  In  the 
following  discussion,  the  77  direction  (or  y  direction  in  the  rectilinear  case)  is  also  referred  to 
as  the  north-south  direction,  while  the  ^  (or  x)  direction  is  associated  with  the  east-west 
direction.  The  user  currently  has  the  following  options  with  regard  to  boundary  conditions: 
C  (or  x)  direction  77  (or  y)  direction 

•  periodic  •  periodic 

•  east  wall/ west  wall  •  open  (north  and  south) 

•  east  wall/west  open 

Pree-slip  no-normal-flow  boundary  conditions  are  applied  at  the  weJls.  Gravity-wave  radi¬ 
ation  is  approximated  at  the  open  boundaries. 

In  general,  boundary  conditions  can  be  treated  in  one  of  two  ways;  (i)  one  can  appro¬ 
priately  assign  values  to  points  “outside”  the  computational  domain  and  apply  the  same 
algorithm  used  in  the  interior,  or  (ii)  one  can  apply  a  different  algorithm  at  the  boundary. 
The  treatment  for  walled  boimdaries  in  a  rotating  flow  follows  the  second  approach,  while 
the  others  use  the  first  approach.  Recall  that  the  r-th  order  ENO  scheme  requires  an  r  -t- 1 
adaptive  stencil,  so  to  follow  the  first  approach  r  -t- 1  points  outside  the  domain  must  be 
assigned. 

In  the  discussion  that  follows,  consider  a  computational  grid  that  is  discretized  into 
M  X  N  grid  points, 

G,  i  =  0,...,M-l 

j  =  0,...,N  -1. 

The  western  and  eastern  boundaries  are  located  at  ^0  and  Cm— i?  respectively.  The  southern 
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and  northern  boundaries  are  located  at  r)o  and  77jv-i,  respectively. 


7.1  Periodic  BCs 

The  implementation  of  periodic  boundary  conditions  is  trivial.  In  the  77  direction,  for 
example,  we  simply  set 

g-j  =  Qn-j 

QN-i+j  =  gj-u  for  ;  =  l,...,r  +  l. 

7.2  Radiation  BCs 

Radiation  boundary  conditions  are  easily  implemented  in  the  model  since  the  ENO  scheme 
is  a  characteristic-based  scheme.  In  the  eigenspace,  the  sign  of  the  eigenvalue  indicates 
the  direction  of  wave  propagation  along  the  characteristic.  When  the  wave  propagation  is 
directed  out  of  the  domain,  the  radiation  treatment  calls  for  extrapolation  of  the  Riemann 
invariant  to  the  grid  points  “outside”  the  domain.  When  the  wave  propagation  is  directed 
into  the  domain,  the  value  of  the  Riemann  invariant  outside  the  domain  is  prescribed.  In 
the  current  model  implementation,  the  prescribed  values  for  the  incoming  waves  are  the 
initial  conditions  along  the  boundary. 

In  the  rj  direction,  for  example,  we  obtain  the  Riemann  invariants  by  projecting  q  using 
the  left  eigenvectors  oi  B  =  dG/dq^  i.e.. 


1 

0 

— n 

Tt  =  p-^g, 

p-i  = 

0 

1 

2c 

1 

1 

lo 

1 

2c 

Tc{v  +  c) 

The  eigenvalues  of  jB, 

( ,  A^^^ ,  A^^^ )  —  {v^v  c^v  —  c) 

reveal  the  direction  of  wave  propagation  at  the  northern  and  southern  boundaries  and 
determine  how  the  value  of  the  Riemann  invariants  outside  the  domain  will  be  set.  At  the 


31 


southern  boundary  (at  grid  point  i  =  0),  for  example,  the  algorithm  for  the  p-th  Riemann 
invariant  is, 

if  >  0  then 

t)  =  (770,  i  =  0),  ;•  =  1, . . . ,  r  +  1 

else 

(7?o,  t),  i  =  1, . . . ,  r  + 1. 

The  value  of  the  state  vector  outside  the  domain  is  then  set  by  projecting  back  to  physical 
space  via  the  right  eigenvectors  of  B,  i.e., 

^  1  It  u  ^ 

g  =  PIZ.,  P  =  0  V  +  c  V  -  c 

\0  1  I  ) 

7.3  Walled  BCs  for  Non-rotating  Flows 

In  the  absence  of  rotation,  u  =  0  implies  =  0  and  =  0  (see  Equation  (7)).  Therefore 
the  no-normal-flow  condition  at  the  western  boimdary  (grid  point  t  =  0),  for  example,  can 
be  satisfied  by  simply  setting, 

tto  =  0,  and 

U—i  =  —Ui 

V-i  =  Vi 

h-i  =  hi,  for  z  =  1, . . .  ,r -t- 1. 

In  this  scenario,  the  points  outside  the  domain  are  typically  referred  to  as  image  points  or 
ghost  points. 

7.4  Walled  BCs  for  Rotating  Flows 

While  the  use  of  ghost  points  works  well  in  the  non-rotating  case,  an  alternative  approach 
is  necessary  when  rotation  is  present.  In  rotating  flows,  u  =  0  implies 

fov=^^F-\.  , 

mi 
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Even  though  we  could  still  anti-image  u  at  the  wall  (e.g.,  u-i  —  —Ui),  we  do  not  know  the 
functional  form  of  h  (or  v),  so  we  cannot  effectively  assign  values  to  h  and  v  at  the  ghost 
points  that  would  maintain  the  geostrophic  relationship  to  a  high  order  of  accuracy.  Ex¬ 
trapolation  from  the  interior  to  the  ghost  points  was  tested  (up  to  4th-order  extrapolation), 
but  a  mismatch  in  the  truncation  error  between  the  ghost  point  and  the  wall  point  resulted 
in  an  error  in  the  flux  at  the  wall  that  eventually  contaminated  the  numerical  solution.  For 
rotating  flows,  we  have  not  found  a  satisfactory  method  to  update  the  flow  field  at  the  wall 
using  ghost  points.  Instead,  we  have  implemented  a  boundary  treatment  that  locally  solves 
a  Riemann  problem  to  update  the  grid  points  at  the  wall  using  only  interior  information. 
Consider  the  model  equations  (see  Equation  (10))  rewritten  as, 

q^  +  ±[F{q)]^  =  Q 

where 

Q  =  -J-[G{q)]r,  +C  +  T  +  T>  +  M. 

m2 

Applying  the  projection  matrix  to  the  system,  i.e.. 


1 


S-^Qt  +  —{S-^AS)S-^qc  =  S-^Q 

mi 


A  = 


dF 


aq’ 

yields  the  characteristic  form. 


S-^  =  I  Z(2)  1  = 
f(3) 


1 


/  0  1 


\ 


^  0  -^(«-c) 

V  0  ^(^■i'^)  J 


Rt  H - A.Rc  —  Q 

mi 


where 


(  ^  \ 

(  u  \ 

R  = 

« -1-  2c  j  ,  A  = 

u  +  c 

^u  —  2cj 

1  u-c  y 
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and 


V  e®-(»+c)c‘®'  7 


(see  Section  4).  Since  u  =  0  at  a  wall,  we  have 


=  QP'^ 

mi  ^ 

iP'^-—rP  = 

mi  ^ 


Note  that  since  when  u  =  0,  the  flow  at  the  wall  can  be  obtained  from  the 

interior  flow  fields,  and  no  information  outside  the  domain  is  necessary.  For  the  eastern 
wall  at  C  =  Cw-i?  we  solve  for  R^^^  and  for  the  western  wall  at  C  =  Co,  we  solve  for 
and  R^^\ 

Consider  the  wall  at  the  eastern  boundary,  C  =  Cm-i,  which  we  now  denote  by  Cw  The 
characteristic  quantities  at  the  wall  at  time  t",  and  are  computed  by  locating 

the  characteristic  that  intersects  the  wall  at  the  required  time  level,  as  depicted  in  Figure  2. 
(The  vector-component  index  has  been  moved  to  the  left  of  the  vector  variable  for  notational 
clarity.)  The  characteristic  quantities  at  the  wall  are  advanced  in  time  in  a  manner  that  is 
consistent  with  the  time-stepping  scheme  used  in  the  interior.  For  example,  when  the  3rd- 
order  TVD  Runge-Kutta  is  used  in  the  interior  (see  Section  6),  the  characteristic  quantities 
at  the  wall  are  advanced  according  to. 


=  ^^'>Rl  +  AP^Ql 
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=  ^^'^ri  +  y{!-''^Qw+^^^qi  +  ^^^^qV}- 

(It  is  easily  verified  that  this  form  of  the  Runge-Kutta  is  equivalent  to  that  defined  by 
Equations  (34)-(36).)  The  computation  of  involves  the  evaluation  of  R{C*)- 

Here  C*  is  the  interior  location  of  the  chaxacteristic  that  intersects  the  wall  at  the  next  partial 
time  step  (see  Figure  2)  and  is  given  by, 


c 

c 

d 


c _ La*c” 

mi 

r _ L^(c»i+c“) 

C _ L^/c"  +  c“  +  Ac^ ) 

mi  D 


Four-point  Lagrange  interpolation  is  then  used  to  evaluate  from  ^'^^Rw,  ^“^^Rw-i, 

^'^^Ru,-2,  and  ^'^'>Ru}-3.  The  new  values  of  R  are  projected  back  using  the  right  eigenvectors 
of  A  =  ^  to  update  the  physical  flow  variables  q. 

Since  the  values  at  the  wall  are  updated  by  this  alternative  method,  the  flux  at  the  half¬ 
grid  point  adjacent  to  the  wall  does  not  have  the  same  truncation  error  as  fluxes  computed  at 
points  farther  in  the  interior  using  the  ENO  algorithm,  and  therefore  high-order  accuracy 
of  the  derivative  is  not  guaranteed.  To  help  alleviate  this  slight  mismatch,  the  first  two 
interior  points  adjacent  to  the  wall  are  weakly  smoothed  after  the  flow  has  been  advanced 
to  the  new  time  level. 


0"+’: 


_  ^71+1  I  c  I  _n4*l 

i  — 1  ^2Qyj—i  i+l 


for  i  =  1,2. 


5l  +  52  +  53 

For  the  original  application  of  the  model,  I  typically  used  the  smoothing  coefficients. 


(S1)«2,S3)  =  (1)38,1). 
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8  Model  Preparation  and  Execution 


In  this  section,  instructions  are  provided  on  how  to  prepare  and  execute  model  runs.  In  the 
discussion  that  follows,  the  reader  may  find  it  helpful  to  refer  to  the  source  code  since  it 
contains  comments  as  well. 

The  model  code  is  written  in  Fortran  77  and  uses  C-directives  (i.e.,  #def  ine  and  #ifdef 
constructs)  to  organize  the  somce  code  corresponding  to  various  model  options.  The  source 
code  is  divided  into  four  files:  two  header  files,  enores.h  and  gpath.h,  a  common  file, 
enoswcom.f ,  and  the  main  program,  enosw.F. 

8.1  Header  File  enores  .h 

Within  the  header  file  enores  .h,  the  user  must  set  the  parameter  specifications  for  the  grid 
size  and  the  numerical  accuracy  of  the  scheme. 

•  Q1  specifies  the  spatial  order  of  accuracy  for  the  ENO  scheme.  The  ENO  scheme  is 
intended  to  be  a  high-order  numerical  scheme.  For  instance,  one  would  not  specify 
Ql=l  (specifying  a  2-point  stencil)  since  this  would  greatly  restrict  the  stencil  selection 
process  which  forms  the  cornerstone  of  the  ENO  algorithm.  The  setting  Ql=3  (4-point 
stencil)  is  more  typical. 

•  Q2  specifies  the  temporal  order  of  accuracy  for  the  Runge-Kutta  scheme.  Runge- 
Kutta  schemes  for  Q2=2  and  Q2=3  are  implemented  (see  Equations  (32)-(33)  and 
Equations  (34)-(36)). 

•  M  and  N  specify  the  grid  size. 

•  IRSIZ=2*M*N  currently  specifies  the  record  length  for  the  direct-access  unformatted 
(binary)  I/O  of  one  double-precision  MxN  model  data  field. 
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8.2  Header  File  gpath.h 


Within  the  header  file  gpath.h,  the  user  must  provide  information  regarding  the  grid  met¬ 
rics.  If  the  grid  is  rectilinear  with  constant  grid  spacing  in  the  x  (=  C)  and  y  (=  r/)  directions, 
then  Ax  and  Ay  need  to  be  provided  by  setting  the  variables  dxO  and  dyO,  respectively, 
e.g., 

dxO  =0.1 
dyO  =0.1 

The  variable  gpath  in  this  case  is  inactive  and  is  used  only  in  the  creation  of  a  log  file  at 
the  end  of  the  model  run  (see  below). 

If  a  more  general  orthogonal  curvilinear  grid  is  to  be  used,  then  the  character  variable 
gpath  specifies  the  path  to  the  grid  metric  data.  The  grid  metric  data  are  expected  in  a 
direct-access  unformatted  (binary)  file  called  swgrid.met,  in  record  1.  The  read  statement 

within  subroutine  Init  of  enosw.F  is, 
c  *  Get  grid  metrics 

open(unit=10 ,f ile=gpath(  1 : Inblnk (gpath) ) // ’ swgrid.met ’ , 

form=’ unformatted’ ,access=’direct’ ,recl=2*IRSIZ,status=’old’) 
read(10,rec=l)  ((hl(i,j),  h2(i,j),  i=0,M-l),  j=0,N-l) 

close (10) 

in  which  the  arrays  hl(i,j)  and  h2(i,j)  store  the  value  of  the  grid  metrics  mi  and  m2 
(see  Equations  (5)-(6)). 

8.3  Common  File  enoswcom.f 

The  file  enoswcom.f  defines  the  common  blocks  and  contains  all  of  the  C-directives  con¬ 
trolling  the  various  model  options,  described  below.  The  user  sets  the  various  C-directives 
by  un-commenting  the  #def  ine  statements  corresponding  to  the  desired  options  and  com¬ 
menting  out  the  others. 

•  ENO  algorithm.  As  discussed  in  Section  3.1,  the  Hybrid  ENO-Roe  algorithm  is  the 
recommended  algorithm.  Therefore,  the  recommended  settings  are 


Sdefine  HYBRID 
#define  ROE 

The  Hybrid  ENO-LLF  algorithm  might  also  be  considered, 

ttdefine  HYBRID 
#defiiie  LLF 

but  in  my  experience,  the  ENO-Roe  solutions  have  been  superior  to  those  produced 
using  ENO-LLF. 

•  Computational  grid  type.  If  the  grid  is  uniform,  un-comment  the  statement, 

#define  UHIGRID 

When  UNIGRID  is  “on”,  the  grid  spacing  must  be  specified  in  the  header  file  gpath.h 
(see  above).  If  the  grid  is  not  uniform,  it  is  expected  that  the  grid-metric  input  file 
can  be  found  in  the  location  indicated  within  gpath.h  (see  above). 

•  Bottom  stress  parameterization.  A  linear  bottom  stress  is  selected  by  turning  on 
LSTRESS.  If  LSTRESS  is  on,  the  linear  drag  coefficient  f  will  be  taJsen  firom  the  pa¬ 
rameter  file  enosw.parms  (see  below).  At  present,  two  nonlinear  drag  coefficients  are 
coded.  If  BSTRESS  is  on,  the  drag  coefficient  derived  by  Large  and  Pond  will  be  used; 
if  CWSTRESS  is  on,  the  formula  derived  from  the  Coastal  Waves  96  data  will  be  used. 
(See  Section  5.) 

•  Forcing.  During  the  initialization  phase  of  the  model,  the  steady  component  of  the 
forcing  field  (see  Section  5)  is  read  from  the  input  fide  enosw_in.dat,  described  below. 
If  the  C-directive  STDYFORC  is  on,  the  forcing  is  assumed  to  be  steady  (i.e.,  the  steady 
component  is  the  total  field) . 
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If  there  is  a  time-dependent  component,  it  should  be  coded  in  subroutine  Forcing  of 
enosw.F  and  the  C-directive  STDYFORC  should  be  turned  off.  In  the  original  applica¬ 
tion  of  the  model,  time-dependent  forcing  was  used  and  the  computational  grid  was 
curvilinear.  The  time-dependent  pressure  gradient  forcing  was  defined  on  the  x-y  grid 
and  then  mapped  to  the  ^-77  grid.  Therefore  in  the  current  version,  if  STDYFORC  is 
on  and  UNIGRID  is  off,  the  program  will  attempt  to  read  the  curvilinear  grid  points 
and  the  (x,y)-to-(C,eta)  transformation  matrix  (see  Section  8.7)  in  subroutine  Init 
of  enosw.F  with  the  statements, 

c  *  Get  the  grid  for  the  HARD-CODED  pressure  forcing  function 
open(unit=10 , f ile=gpath( 1 :  Inblnk (gpath) ) // ’ swgrid . pts ’ , 

.  form=’unformatted’ ,status=’old’ ) 

read(lO)  ((x(i,j),  y(i,j),  i=0,M-l),  j=0,N-l) 

close(lO) 

c  *  Get  the  (x,y)->(zeta,eta)  mapping  matrix  to  map  grad(P) 
open (\mit=10 , f ile=gpath ( 1 : Inblnk ( gpath) ) // ’ swgrid . map ’ , 

form=’ unformatted’ ,access=’ direct’ ,recl=4*IRSIZ,status=’old’) 
read(10,rec=2)  ((rotillU,  j)  ,rotil2(i,  j)  , 

roti21(i,j) ,roti22(i, j) ,i=0,M-l) ,  j=0,N-l) 

close(lO) 


•  Boundary  conditions.  A  number  of  boundary  conditions  have  been  implemented. 
Refer  to  Section  7  for  a  complete  description.  For  example,  to  specify  a  periodic 
channel  in  a  rotating  system,  the  C-directives  should  be  turned  on  and  off  as, 

c#define  XPERIODIC 
c#define  East GHOST 
c#define  WestGHOST 
#define  EastWALL 
#defiiie  WestWALL 
c#define  WestCHAR 
#define  YPERIODIC 
c#define  YCHAR 
c#define  YORLANSKI 


For  a  channel  in  a  non-rotating  system  with  radiation  conditions  at  the  channel  ends, 
the  settings  would  be, 


40 


c#define  XPERIODIC 
#define  EastGHOST 
#define  WestGHOST 
c#define  East WALL 
c#define  WestWALL 
c#define  WestCHAR 
c#define  YPERIQDIC 
#define  YCHAR 
c#define  YORLANSKI 


Note  that  the  Orlanski  radiation  treatment  is  implemented  for  comparative  purposes 
only.  The  YCHAR  option  should  be  used  for  radiation  boundary  conditions. 

•  Smoothing  for  East  WALL  or  WestWALL. 

If  EastWALL  or  WestWALL  is  on,  the  C-directive  SMOOTHI  should  also  be  on  (see 
Section  7).  The  smoothing  coefficients  are  provided  in  the  input  parameter  file 
enosw.parms,  described  below. 

8.4  InitialConditionDataFileenosw_in.dat 

It  is  assumed  that  the  initial  flow  fields  (i.e.,  the  flow  velocities,  u  and  v,  and  the  layer 
thickness  h)  and  the  steady  component  of  the  pressure-gradient  forcing  (i.e.,  VP°)  reside  in 
a  direct-access  unformatted  file  called  enosw_in.dat.  The  read  statement  within  subroutine 
Init  of  enosw.F  is, 

c  *  Get  initial  u,v,h,  and  forcing  fields 

open(unit=10,file=’enosw_in.dat’  ,fonn=’unformatted’ , 
access= ’ direct ’ ,recl=5*IRSIZ,status=’old’) 
read(lO)  ((u(i,j ,0) ,v(i, j ,0) ,h(i,j ,0) ,pxfO(i, j) ,pyfO(i,j) , 

i=0,M-l),  j=0,N-l) 

in  which  the  following  correspondence  is  made, 

u(i,j,0)  -H  «(Ci,?7j,<  =  0) 
v(i,j,0)  4-  v{Ci,r]j,t  =  0) 

h(i,j,0)  ■(-  h{Ci,r]j,t  =  0) 
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pxfO(i,j)  ■<— 

pyfO(i,j) 


1  dP^iCv) 
mi  dC 
1  dP^{C,v) 
m2  dr] 


After  reading  the  initial  flow  fields,  the  variables  u(, ,)  and  v(, ,)  are  reassigned  to  hold 
the  flux  variables  U  -uh  and  V  =  vh,  respectively. 

8.5  Model  Parameter  File  enosw.parms 

A  run-time  parameter  file,  enosw.parms,  is  also  supplied  to  the  model.  An  example  of  this 
ASCII  file  is  shown  below. 


10.  10.  500.. 

0.1  1.0 
0.0 
0.01 

0.  200.  50. 

1.  38.  1. 


L*  (km) ,  U*  (m/s) ,  D*  (m) 
f 0=l/Rossby=fL/U,  Fr-2=1/ (Froude~2) 
r_linear  (ignored  if  BSTRESS/CWSTRESS  is  on) 
dt 

t start,  tfinal,  tdump 
smoothing  coefficients  (sl,s2,s3) 


•  Line  1  specifies  the  length,  velocity  and  depth  scales  that  will  be  used  in  the  hard¬ 
coded  formulas  that  parameterize  the  nonlinear  bottom  stress  in  subroutine  Friction 
when  either  BSTRESS  or  CWSTRESS  is  on  (see  Section  5). 


•  Line  2  specifies  the  inverse  Rossby  number,  /o  =  fL*/U*,  and  scaling  inverse-squared 
Proude  number,  =  9'D*/U*^  (see  Section  2). 

•  Line  3  specifies  the  hnear  (nondimensional)  drag  coefficient  if  the  C-directive  LSTRESS 
is  on.  If  a  nonlinear  bulk  formula  is  specified  (i.e.,  C-directive  BSTRESS  or  CWSTRESS 
is  on),  then  the  velocity-dependent  value  for  the  drag  coefiicient  is  hard-coded  in 
subroutine  Friction  and  the  value  of  the  linear  drag  coefiicient  in  enosw.parms  is 
ignored  (see  Section  5). 

•  Line  4  specifies  the  nondimensional  time  step.  At  (see  Section  6). 
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•  Line  5  specifies  the  time  associated  with  the  beginning  of  the  simulation  (tstart), 
the  time  associated  with  the  end  of  the  simulation  (tf  inal),  and  the  time  interval  at 
which  the  model  data  will  be  output  (tdump). 

•  Line  6  specifies  the  smoothing  coefficients  that  are  required  when  characteristic-based 
walled  boundary  conditions  are  in  effect  (C-directive  EastWALL  and/or  WestWALL)  (see 
Section  7). 

8.6  Execution 

To  conduct  model  runs, 

1.  Select/generate  a  computational  grid  and  edit  gpath.h  appropriately. 

2.  Set  the  parameters  in  enores.h. 

3.  Set  the  C-directives  that  control  the  model  options  in  enosw.parms. 

4.  Generate  enosw_in.dat,  the  initial-condition  data  file. 

5.  Edit  enosw.parms  as  necessary. 

6.  Make  enosw. 

7.  Execute  enosw.  (There  are  no  command-line  arguments.) 

8.7  Output  files 

The  executable,  enosw,  generates  the  output  files,  u.dat,  v.dat,  and  h.dat,  containing  the 
data  for  u,  v,  and  h,  respectively.  If  the  C-directive  STDYFORC  is  off,  two  additional  out¬ 
put  files,  px.dat  and  py.dat,  are  created  for  the  time-dependent  pressure-gradient  forcing 
field.  Each  file  is  a  direct-access  unformatted  file.  For  example,  u .  dat  is  opened  with  the 
statement, 

open (unit=20,file=' u.dat’  ,form=’mformatted’ , 
access= ’ direct ’ ,recl=IRSIZ) 
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The  initial  data  are  placed  in  record  1,  and  the  remaining  records  correspond  to  time  tdump, 
time  2*tdump,  time  3*tdump,  etc.,  where  tdump  is  specified  in  enosw.parms.  For  example, 
the  u-velocity  at  time  t  =  0  is  output  with  the  statement, 

write (20, rec=l)  ((u(i,j ,0)/h(i, j ,0) ,  i=0,M-l),  j=0,N-l) 

(Recall  that  the  variable  u( , ,)  holds  the  fiux  U  =  uh.) 

A  log  file,  enosw.log,  is  also  created  to  catalogue  the  input  parameters  that  were 
specified  for  the  model  run.  An  example  is, 

(M,N)=(100,100)  (Q1,Q2)=(3,3)  dt=0. 01000 

(to ,  tmax ,  tdun^))  =  (  0 . 000 , 200 . 000 , 50 . 000) 

L*(km)=  10.  U*(m/s)=  10.  D*(m)=500. 

f0=  0.100  r=0.0000  Fi=  1.000 

Gridruniform 

(dx,dy)=(  0.100,  0.100) 

Smoothing  coefficients:  1.  38.  1. 

ROE-EF  HYBRID  EastWALL  WestWALL  YCHAR  LSTRESS  STDYFORC  SMOOTHI 
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Appendix:  Curvilinear  Grid-Generation  Program  swgrid.f 


The  shallow-water  model  is  designed  to  approximate  a  solution  on  a  user-specified  orthog¬ 
onal  curvihnear  grid.  If  the  desired  grid  is  rectiUnear  with  uniformly-spaced  grid  points  in 
the  X  and  y  direction,  then  the  user  need  only  specify  Arc  and  Ay  in  addition  to  the  number 
of  grid  points  in  each  direction.  If  a  more  general  orthogonal  grid  is  to  be  used,  the  user 
must  provide  the  (spatially-varying)  grid  metrics  to  the  shallow-water  model. 

One  orthogonal  curvihnear  grid  generation  program,  swgrid.f,  is  included  in  the  model 
package  as  an  example.  In  this  section,  we  review  the  basics  of  coordinate  transformations 
and  provide  a  brief  description  of  swgrid.f  and  how  it  interfaces  with  the  shallow- water 
model. 

If  (C,  y)  are  the  coordinates  in  the  orthogonal  curvilinear  system,  then  the  change  in  the 
position  vector  x  =  {x,y)  =  in  the  Cartesian  system  can  be  written  as, 

6x  =  mi6(^  C  +  “^267}  f) 


where  mi  and  m2  are  the  coordinate  metrics  given  by 


A  vector  {u,v)  on  the  (x,y)-grid  can  be  transformed  to  the  (C,y)-grid  via. 


X^/imiA)  Xr,/im2A)  \ 
^  Y^/imiB)  Yr,l{m2B)  ) 


(u\ 

[v  J 


where 


A  = 

B  = 


(38) 

(39) 


(40) 
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Similarly,  a  vector  {u,v)  on  the  {(^,ri)-grid  can  be  transformed  to  the  («,y)-grid  via, 

j  /  YjjTniA  —XrjmiB 

X^m2B 

The  metrics  m;  and  m2  (Equations  (38)-(39))  must  be  provided  to  the  shallow- water 
model.  The  coordinate  transformation  (40)  can  be  used  during  the  model  initialization  to 
map  flow  fields  specified  on  a  rectilinear  (e.g.,  north-south,  east-west)  coordinate  system 
to  the  model’s  curvilinear  coordinate  system.  During  the  post-processing  phase,  the  inverse 
mapping  (41)  can  be  applied  to  visualize  the  model  output  in  the  rectilinear  coordinate 
system  if  so  desired. 

The  grid-generation  program  swgrid .  f  was  originally  developed  by  Wilkin  and  modified 
successively  by  R.  Signell,  by  R.  Samelson  and  by  A.  Rogerson.  In  short,  the  user  specifies 
the  desired  boundary  and  provides  an  initial  distribution  of  grid  points  along  the  boundary. 
The  orthogonal  curvilinear  grid  is  obtained  by  iteratively  applying  a  conformal  mapping 
algorithm  to  the  gridded  domain.  Detailed  aspects  of  the  algorithmic  approach  and  the 
implementation  will  not  be  discussed  here.  Rather,  a  brief  description  of  how  to  use  this 
program  is  provided  below. 

Several  parameters  must  be  set  in  swgrid. f.  The  first  two  specify  the  grid  size  and  are 
set  by  the  parameter  statement 

parameter  (L=100,M=200) 

located  at  the  head  of  program  swgrid  and  the  subroutines  spline,  cof  x,  and  cof  y.  The 
third  is  the  number  of  iterations  to  be  performed  to  obtain  the  grid,  set  by  the  data  state¬ 
ment 

data  itmax  /15/ 

in  the  first  portion  of  swgrid.  Following  immediately  after  the  data  statement  for  itmax  is 
the  data  statement 
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data  kb2/4/ 


This  parameter  specifies  which  of  the  four  boundaries  will  maintain  its  original  distribution 
of  grid  points  after  the  mapping. 

The  user  specifies  the  geometric  shape  and  the  initial  distribution  of  grid  points  for  the 
western,  southern,  eastern,  and  northern  boundaries  within  subroutine  zl,  and  entry  points 
z2,  z3,  and  z4,  respectively.  The  boundary  position  data  are  stored  in  a  single  complex¬ 
valued  vector  variable  z(i)=(x(i)  ,y(i)),  which  holds  the  position  of  the  grid  point  at 
the  northwest  corner  of  the  domain  in  z(l)  and  the  remaining  data  points  in  successive 
storage  locations,  proceeding  counter-clockwise  aroimd  the  boundary.  In  subprograms  zl, 
z2,  z3,  and  z4,  the  boundary  data  are  defined  parametrically  through  the  variable  s  which 
varies  from  zero  to  one  along  each  portion  of  the  boundary.  In  the  current  setup,  the 
domain  is  30  units  long  in  the  x  direction  and  50  units  long  in  the  y  direction  (variables 
XL  and  YL  in  subroutine  zl)  and  is  discretized  into  L  x  M  =  60  x  100  grid  points.  The 
eastern  boundary  consists  of  a  series  of  bends;  phi(j)  are  the  bend  angles,  measured  from 
due  south,  ybend(j)  are  the  y  positions  of  the  bends,  rc(j)  are  the  radii  of  curvature 
for  the  bend.  The  formula  for  the  eastern  boundary  points  is  obtained  after  a  little  bit  of 
trigonometry.  The  western,  northern,  and  southern  boundaries  are  straight.  Grid  points 
along  the  southern  boundary  are  equally-spaced  while  those  along  the  western  and  northern 
boundaries  are  clustered  non-uniformly.  The  clustering  in  the  current  implementation  is 
achieved  by  mapping  the  parametric  variable  s  to  a  piecewise  continuous  cubic  polynomial. 
Additional  details  are  provided  by  the  comments  within  the  source  code. 

When  the  deformation  of  the  boundary  is  severe,  the  grid  that  is  generated  may  not  be 
orthogonal  near  the  edges  of  the  domain.  Two  additional  parameters  have  been  introduced 
at  the  beginning  of  program  swgrid  to  chp  the  grid  at  the  northern  and/or  southern  ends 
during  output. 
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parameter  (jclipn=0,  jclips=5) 


In  this  case,  the  grid  points  indexed  by  j  <  5  are  eliminated. 

Program  swgrid.f  produces  five  output  files: 

•  swgrid.met,  a  direct-access  unformatted  (binary)  file  containing  the  grid  metrics,  m\ 
and  7712,  on  the  Sadourney  C  grids.  The  metrics  on  the  /i,  ti,  and  v  grids  axe  stored  in 
records  1,  2,  and  3,  respectively.  The  shallow- water  model  accesses  the  grid  metrics 
on  the  h  grid  (record  1). 

•  swgrid.map,  a  direct-access  unformatted  file  containing  the  coordinate  transformation 
metrics  that  appear  in  Equations  (40)  and  (41).  The  (a;,y)-to-(C,77)  transformation 
matrix  is  stored  in  record  1.  The  (C,7?)-to-(x,y)  transformation  matrix  is  stored  in 
record  2. 

•  swgrid.pts,  an  unformatted  file  containing  the  grid  points  (x,  y)  =  (X(C,  rj),  Y (C,  »?))• 

•  swgrid.bdry,  an  ASCII  file  containing  the  boundary  points. 

•  mesh.dat,  an  ASCII  file  to  ease  the  plotting  of  the  grid  mesh. 

Note  that  for  the  direct-access  binary  files,  the  record  length  is  the  minimum  required  for 
double-precision  output. 
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